First-Order Phase Transition with Breaking of Lattice Rotation Symmetry in 
Continuous-Spin Model on Triangular Lattice 
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Using a Monte Carlo method, we study the finite-temperature phase transition in the two- 
dimensional classical Heisenberg model on a triangular lattice with or without easy-plane 
anisotropy. The model takes account of competing interactions: a ferromagnetic nearest- 
neighbor interaction Ji and an antiferromagnetic third nearest-neighbor interaction J3. As 
a result, the ground state is a spiral spin configuration for —4 < J1/J3 < 0. In this structure, 
global spin rotation cannot compensate for the effect of 120-degree lattice rotation, in contrast 
to the conventional 120-degree structure of the nearest-neighbor interaction model. We find 
that this model exhibits a first-order phase transition with breaking of the lattice rotation 
symmetry at a finite temperature. The transition is characterized as a Z2 vortex dissociation in 
the isotropic case, whereas it can be viewed as a Z vortex dissociation in the anisotropic case. 
Remarkably, the latter is continuously connected to the former as the magnitude of anisotropy 
decreases, in contrast to the recent work by Misawa and Motome [J. Phys. Soc. Jpn. 79 (2010) 
073001.1 in which both the transitions were found to be continuous. 
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1. Introduction 

In frustrated magnets, there is no ground state wiiere 
all interactions are satisfied energetically. Consequently, 
conventional magnetic orders are suppressed. ^^^-^ The 
frustration appears in antiferromagnets on triangle-based 
lattices such as a triangular lattice, a kagome lattice, and 
a pyrochlore lattice. It is well-known that interesting or- 
dering phenomena appear in frustrated magnets. For ex- 
ample, order by disorder on a triangular lattice,*'^' as 
well as on a kagome lattice, has been reported. In par- 
ticular, antiferromagnets that have a triangular lattice 
are important, because they represent simplest realiza- 
tion of geometrical frustration. Thus, the antiferromag- 
netic model on a triangular lattice has been studied ex- 
tensively. ^"^-'^^ 

The chalcogenide insulator NiGa2S4 has recently been 
investigated by Nakatsuji and his collaborators. ^^"^^^ 
NiGa2S4 has an almost perfectly equilateral triangular 
lattice that consists of Ni^+ ions {S = 1). From the sus- 
ceptibility result, the Weiss temperature was estimated 
to be 9vi = —80 K, which means that strong antiferro- 
magnetic (AF) correlation exists. This strong AF cor- 
relation caused by a third nearest-neighbor (NN) in- 
teraction has been suggested from photoemission spec- 
troscopy measurements,^^) neutron-scattering measure- 
ments,^^'' and first-principle calculations.^^) According to 
these studies, the second NN interaction, which is much 
weaker than the first and third NN interactions, is neg- 
ligible. The magnetic specific heat exhibits an double- 
peak structure centered near the Weiss temperature |0w| 
and 10 K, and the peaks do not diverge. No mag- 
netic long-range order was observed down to 0.35 K in 
experiments. Meanwhile, short-range incommensurate 
(IC) correlation, the correlation length of which is a few 
sites, was detected below 10 K. Furthermore, a linearly- 



dispersive gapless coherent mode in two dimension was 
observed in this material. Thus, the magnetic short- 
range order is expected to be strongly related to two- 
dimensional properties. 

Initially, NiGa2S4 was considered to be an AF material 
on a triangular (AFT) lattice, but it has been found 
that simple AFT Heisenberg models^^'^^^ do not describe 
NiGa2S4 well. One discrepancy is that the correlation 
length of spins monotonically develops in the simple AFT 
models, whereas it is limited within the range of 10 lattice 
constants in NiGa2S4. Another discrepancy is that the 
magnetic short-range ordered phase does not exist in the 
simple AFT models. 

To explain these discrepancies, several theoretical 
models and scenarios have been proposed, (i) The S = 
1 quantum AF Heisenberg model with NN bilinear- 
biquadratic interactions on a triangular lattice^'^"^^' de- 
scribes the short-range magnetic order, and the model 
can explain the gapless linear dispersion observed in 
NiGa2S4. Although the quantum bilinear-biquadratic 
model succeeded in describing a couple aspects of 
NiGa2S4, the development of nematic order in the model 
is not consistent with experimental observations, (ii) The 
Z2 vortex scenario has been proposed to describe the 
phase transition where the specific heat does not diverge. 
In the Z2 vortex scenario, a Z2 topological vortex ex- 
plains the topological phase transition in the classical 
Heisenberg model with AF NN bilinear interaction on 
a triangular lattice (AF Ji model).*' However, an 
IC correlation has been observed in NiGa2S4, and the AF 
Ji model shows a commensurate 120-degree structure. 

The AFT model is known to exhibit an IC phase when 
the model has an AF NN interaction and an AF second 
NN interaction.^^' Therefore, since an IC phase ap- 
pears through competition among several types of inter- 
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actions, the classical Heisenbcrg model on a triangular 
lattice with NN and third NN bilinear interactions ( Ji- 
J3 model) is expected to be a promising model for de- 
scribing NiGa2S4. Recently, we studied the Ji- J3 model 
with the interaction ratio J1/J3 = —1/3^°^ and found the 
development of the incommensurate correlation at low 
temperatures. We clarified the occurrence of a, first-order 
phase transition with breaking of the lattice rotation 
symmetry. However, this was inconsistent with experi- 
mental results for NiGa2S4, in which a first-order phase 
transition has not yet been observed. However, the Ji- J3 
model reported in our previous letter suggested an inter- 
esting phcinomcina: a first-order phase transition caused 
by discretization of rotation symmetry. In the previous 
paper, we considered the mechanism of the first-order 
phase transition, but several points remain to be clari- 
fied. 

The aim of this paper is to explain the issues listed 
below, (i) In the previous study, a finite-size scaling anal- 
ysis did not work well because of a correction due to the 
incommensurability. Thus, in this paper, we adjust the 
interaction ratio Ji/ J3 to realize a commcnsmate ground 
state. This adjustment enables us to eliminate an irrele- 
vant but quite large correction. In addition, owing to this 
adjustment, we can check whether the incommensurabil- 
ity is essential for the occurrence of the first-order phase 
transition, (ii) In the Ji- J3 model, there is the Z2 point 
defect. Thus, we confirm the occurrence of the dissocia- 
tion of Z2 vortex pairs in the Ji- J3 model. In addition, 
we investigate whether the first-order phase transition 
and the dissociation of vortex pairs separately occur 
at finite temperatures, (iii) In NiGa2S4, the existence of 
easy-plane anisotropy has been observed through elec- 
tron spin resonance experiments.^^' Additionally, in 
the AF ,/i model with easy-plane anisotropy, different 
types of phase transitions from the isotropic case have 
been observed by Capriotti et al?'^'^^^ and Misawa and 
Motome.^^^ Thus, it is interesting to investigate whether 
the same alteration occurs in the J1-J3 model. Accord- 
ingly, wc! examine whether easy-plane anisotropy changes 
the phase transition in the Ji- J3 model. 

In §2 we introduce the J\-J% model and investigate its 
ground state. We find that the ground state has three 
types of spiral configurations for —4 < J1/J3 < 0. In 
§3 we show the results for the thermal properties in 
the isotropic, Ji-J:\ model. To clarify the irrelevance of 
incommensurability in the first-order phase transition, 
we use the interaction ratio to realize a commensurate 
ground state. It is clear that incommensurability is not 
essential for the first-order phase transition, because the 
first-order phase transition always takes place. Further- 
more, to confirm the occurrence of the dissociation of Zi 
vortex pairs, we calculate the temperature dependence 
of the number density of Z^, vortices. We find that the 
dissociation of Z-2. vortex pairs occurs at the first-order 
phase transition temperature. In §4 we discuss the ef- 
fects of easy-plane anisotropy. We show the results for 
the thermal properties of the XY spin limit in the J1-J3 
model. In this case, we also confirm the existence of the 
first-order phase transition accompanying the Z vortex 
dissociation. Prom easy-plane anisotropy dependence of 



the transition temperature and the latent heat, we find 
that the Z vortex dissociation is continuously connected 
to the Z2 vortex dissociation in the J1-J3 model. In §5 we 
investigate the interaction ratio dependence of the tran- 
sition temperature and the latent heat. From the results, 
we summarize the phase diagram of temperature versus 
interaction ratio and easy-plane anisotropy. Section 6 is 
devoted to discussion and summary. 

2. Model and Ground State 

The Hamiltonian of the J1-J3 model is given by 

'H = Ji Si - Sj + J3 ^ Si - Sj, (2.1) 

where Si is the vector spin of unit length. The first sum 
is taken over NN pairs of sites, and the second sum is 
taken over third NN pairs (see Fig. 1). Here, we adopt 
the periodic boundary condition. 

In this section, we derive the ground state spin con- 
figuration by mean-field calculation. The ground state 
of a classical Heisenberg model is described by a spiral 
configuration with a wave vector k that minimizes the 
Fourier transform of interactions J(fe).^^'^^) The spiral 
configuration is given by 



Sj = ilcos(fc • ri) — /sin(fe • rj, 



(2.2) 



where R and I are two arbitrary orthogonal unit vec- 
tors, and ri is the position of site i in the real space 
on a triangular lattice. Hereinafter, the lattice constant 
is set to unity. We define the primitive translation vec- 
tors of a triangular lattice and its reciprocal lattice as 
ai = (1,0), 02 = (l/2,V3/2), 61 = 27r(l,-l/V3), and 
b2 = 27r(0, 2/V3), respectively. If the NN and third NN 
interactions exist, the Fourier transform of interactions 
J{k) is given by 

J{k)/Js = Ji/J3 I cos(fcx) -h 2cos cos [ -^^1/ ) [ 



+ 



|cos {2kx) -\- 2 cos {kx) cos {V^ky^ | . 



(2.3) 



In this paper, we consider the case where the third NN 
interaction is antiferromagnetic (J3 > 0). When J3 is 
positive, the ground state can be classified into four types 
depending on the interaction ratio J1/J3. 

(i) Ferromagnetic ground state (J1/J3 < —4) 

In this region, the ferromagnetic interaction Ji is 
dominant, and thus J(fe) has a minimum value at 
fc = in the first Brillouin zone. The point at fe = 
is indicated by point 'A' in Fig. 2. The order param- 
eter space of the ferromagnetic state is isomorphic to 
the two-dimensional sphere S2 which corresponds to 
S0(3)/U(l) symmetry. The order parameter space 
is equivalent to the degrees of freedom of the order 
parameter that characterizes the ground state. 

(ii) Spiral ground state (—4 < J1/J3 < 0) 

In this region, neighboring spins in the ground state 
rotate with pitch governed by the minimum point 
of J{k). The shift of the minimum point of J{k) 



3 



from fc = is caused by tlie competition between 
Ji and J3. Unlike tlie AF Ji model, the pitch of the 
spin rotation is not equal to 120 degrees. At the six 
points listed below, J{k) takes its minimum value: 

k = ±(fc,0),±(fc/2, V3fc/2),±(fc/2,-V3fc/2). 

(2.4) 

The wave number k ~ \k\ is given by 

2(sinfc + sin2fc) 

Jihh r-r- — : — n~, (,2-oj 

sm k + sm -^k 

where the range of the value of fc is < fc < 27r/3. A 
fc-point on the x-axis that describes a spiral state is 
indicated by point 'B' in Fig. 2. The point 'B' runs 
from point 'A' to 'C as the value of J1/J3 increases. 
In the case of isotropic Heisenberg spins, the order 
parameter space of this spiral state is isomorphic to 
the three-dimensional real projective space P3 which 
corresponds to S0(3) symmetry.®^ In this case, the 
spin configurations characterized by k and —k can 
be regarded as the same structure. This is because 
the difference between the two spin configurations 
can be eliminated by global SO (3) spin rotation. 
Thus, there are three types of structures, which can 
be characterized by A; = (fc, 0), {k/2,^k/2), and 
(fc/2, — \/3fc/2) at the ground state in the isotropic 
Heisenberg model. These three states are separated 
by energy barriers. In the spiral spin configuration 
characterized by these wave vectors, there are two 
types of rotational pitches. Along one of the three 
axes, spins rotate with the wave number k (axis 1 
in Fig. 1), while along the other axes, spins rotate 
with the wave number k/2 (axis 2 and axis 3 in 
Fig. 1). The three types of structures correspond 
to an axis, characterized by k, being selected from 
among the three axes. In this structure, global spin 
rotation cannot compensate for the effect of 120- 
degree lattice rotation. In the spiral ground state 
region, pitches of the spin rotation can be varied by 
changing the interaction ratio J1/J3. In particular, 
a commensurate spiral spin configuration is realized 
when we set J1/J3 such that k/n is a rational num- 
ber. 

(iii) Four-sublattice 120-degree structure (J1/J3 = 0) 

In this parameter case, the NN interaction Ji is zero, 
and thus the lattice is divided into four indepen- 
dent sublattices. Each sublattice is equivalent to the 
AF Ji model.®' Thus, the ground state is the 
four-sublattice 120-degree structure. The minimum 
points of J(fc) are k — 2tt/3 and k = 47r/3, and these 
points on the x-axis are denoted by points 'C and 
'D' in Fig. 2. The spin configurations characterized 
by eq. (2.4) can be regarded as the same structure. 
This is because these structures are not separated by 
energy barriers in contrast to the case of (ii). The 
order parameter space of the four-sublattice 120 de- 
gree structure is isomorphic to P3. 

(iv) 120-degree structure (J1/J3 > 0) 

In this region, the NN interaction Ji is antiferromag- 
netic, and thus the ground state has the 120-degree 




axis 1 



Fig. 1. (Color online) Schematic of spin configuration at fc = 
(fc,0) and k = 5it/9. The spins of the same color or type of 
arrow form a triangular lattice having twice the lattice constant 
and indicate four sublattices. 

structure. The minimum point of J{k) is located at 
k = 47r/3, and the point on the x-axis is denoted 
by point 'D' in Fig. 2. The order parameter space of 
the 120-degree structure is isomorphic to P3. 

The ground states in the J1-J3 model are summarized in 
Table I. According to this classification, the spiral ground 
state, which has three types of structures, only appears 
in the case of (ii) . Furthermore, the ferromagnetic phase 
and 120-degree structure appear in other models and 
have been thoroughly studied.®''^'''' Thus, the interesting 
parameter range of the J1-J3 model is —4 < J1/J3 < 0. 
In this paper, we investigate the case where the NN inter- 
action Ji is ferromagnetic and |Ji| is smaller than jJsj. 
In §3 and §4, we fix the wave number at fc = Bir/Q. This 
value corresponds to J1/J3 = —0.73425685 = rs^/g from 
eq. (2.5). The ground state at this parameter is the com- 
mensurate configuration. Along one of the three axes, 
the angle between NN spin pairs is 100 degrees (axis 1 
in Fig. 1), while along the other axes, the angle is 50 
degrees (axis 2 and axis 3 in Fig. 1). 

3. Finite Temperature Properties of Heisenberg 
Spin Model 

3.1 First- Order Phase Transition 

To discuss a finite-temperature phase transition, we 
calculate the internal energy per site E and the specific 
heat C by Monte Carlo simulation based on the standard 
heat-bath method. The specific heat at the temperature 
T is defined by 

c = .^M>!, (3,, 

where (• • • ) indicates the thermal average. Hereinafter, 
the Boltzmann constant /cs is set to unity. The interac- 
tion ratio J1/J3 is fixed at ^STr/g in order that the period 
of the ground state becomes eighteen sites. Each simula- 
tion contains 10^-10® Monte Carlo steps per spin at each 
temperature. We conduct 10-60 independent simulations 
for each size to evaluate the statistical errors. We use the 
conditions below to judge convergence. 

(i) Physical quantities calculated by short MC steps 
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Table I. Classification of ground states in the isotropic case. 



interaction ratio 


ground state 


number of structure types 


type of point defect 


J1/J3 < -4 


ferromagnetic state 


one 





-4 < J1/J3 < 


spiral state 


three 


Zi 


Ji/Jz = O(sublattices) 


120-degree structure 


one 


Zi 


J1/J3 > 


120-degree structure 


one 





ky 
k. 



27r/V3 




Fig. 2. Wave vector k that minimizes J(fc) in the first Brillouin 
zone. Bold hexagon indicates the first Brillouin zone, and its 
vertices indicate the wave number of the 120-degree structure. 
Solid lines between vertices of the first Brillouin zone and cen- 
ter point indicate the lines of fc = ±(A;,0), ±(fe/2, \/3fe/2), and 
±(fe/2, —\/Zk/2), respectively. The vertices of dashed hexagon 
indicate the wave number of the 60-degree structure. 

agree with those calculated by long MC steps. 

(ii) Physical quantities observed in simulations started 

from random configurations agree with those in sim- 
ulations started from the ground state. 

The temperature dependence of the internal energy and 
the specific heat are shown in Fig. 3. The specific heat 
exhibits a single peak which becomes sharp as the lat- 
tice size increases. This behavior indicates the existence 
of a phase transition in the Ji- J3 model. The peak po- 
sition moves to lower temperature as the lattice size in- 
creases. The internal energy shows a discontinuous jump 
around the specific-heat peak. This behavior is equiva- 
lent to that of a first-order phase transition observed at 
Ji/J^i ~ —1/3 where the ground state is the IC spiral 
configuration. 

We investigate the energy distribution to confirm 
whether the phase transition is of the first order. Fig- 
ure 4 shows the energy distribution P{E) near the tran- 
sition temperature. A conspicuous bimodal distribution 
appears near the transition temperature reflecting the 
same distribution probability of both paramagnetic and 
spiral states. We conclude that the phase transition is 
of the first order. Even though the ground state is com- 
mensurate, the Ji- J3 model exhibits the first-order phase 
transition. Thus, the IC magnetic structure is irrelevant 
to the first-order phase transition in the Ji- J3 model. 



To estimate the transition temperature Tc for an infi- 
nite system, we adopt the finite-size scaling of the first- 
order phase transition. The finite-size efii'ect of Tc is given 

by 

T,{L) - Tc oc L-"^. (3.2) 

where d is a dimension of the lattice and Tc{L) is the 
transition temperature at the lattice size We esti- 

mate Tc{L) in three ways: (i) the ratio of the existence 
probability of the paramagnetic state P+ and the spiral 
state P- is 1:1, (ii) the ratio P+/P- is 1/2, and (iii) the 
ratio P+/P- is 1/3. The boundary of P+ and P_ is lo- 
cated at the minimum value of P{E) in between the two 
peaks. We calculate P+ and P_ by using the reweight- 
ing method.^^^ Reweighting is carried out by using the 
following relation: 

P(r'; E) oc e-(i/T'-i/T)L<'i5p(j.. (3 3) 

where P(r; E) is the energy distribution at temperature 
T. Figure 5 shows the plot of Tc{L)/J^ versus The 
lines in Fig. 5 are fitting lines calculated by the least- 
squares method. The y-intercepts of the lines are almost 
the same. Therefore, we obtain the transition tempera- 
ture at the limit as L — >■ cxd, 

T,/ J3 = 0.4746(1). (3.4) 

In §5, we discuss the relation between the interaction 
ratio and the transition temperature. 

The evidence of the first-order phase transition can 
be confirmed by other analyses. We define Pmax(-E') and 
Pmin(£') as the lower-energy peak and the minimum of 
P{E) between the peaks at Tc{L), respectively. If the 
first-order phase transition occurs, the finite-size scaling 
of PmaAE) and Pnun{E) is given by-«-40) 

The size dependence of AP is shown in Fig. 6. The al- 
most linear increase of AP with increasing the lattice 
size indicates the existence of latent heat at the limit as 

L — >■ 00. 

The above discussion about the energy distribution 
has been shown that the system has the first-order 
phase transition at a finite temperature. However, the 
development of the correlation length might be nontriv- 
ial. Because the spin configuration in the spiral state 
is anisotropic, the correlation length would exhibit an 
anisotropic development. When an anisotropy is strong 
enough, the critical exponent u of correlation length de- 
pends on the direction. In such a case, it is known that 
a usual finite-size scaling (FSS) analysis is inapplica- 
ble, but a FSS analysis with careful consideration of 
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the anisotropy is required. Since our model pos- 
sibly requires such a special FSS analysis, we check 
the anisotropy of the system by calculating correlation 
lengths which are dependent on the anisotropy. Using 
the Ornstein-Zernike formula,^^^ correlation lengths are 
estimated by ratio of structure factor amplitudes. The 
structure factor S{q) is given by 



Sj)e 



iq-{r 



,) 



and the correlation length is cstiinatc^d l)y 



1 



1, 



(3.6) 



(3.7) 



where <7q is one of six fe's which satisfy cq. (2.4). When 
we assume Qq = (fc, 0), the appropriate q's in eq. (3.7) 
are the three vectors of the following: q-^ = h\/L, 
^2 = 90 ± b2/L, and q-^ = ± (^i + ^2)/-^- Figure 7(a), 
we plot the three correlation lengths $(qi), ^(92)) 
^(^3) for L = 72. Figure 7(a) clearly shows ^(^2) is 
smaller than i,{qi) and ^(93). This indicates that fluctu- 
ations perpendicular to the axis 1 are larger than those 
of the perpendicular to axes 2 and 3 (see Fig. 1). In 
order to examine whether the difference in the fluctua- 
tions brings about a significant effect on the critical ex- 
ponent v, we compare ^(^2) with ^(qfj)(= ^{q^)). The re- 
sult shows that ^(92) differs only a constant factor from 
{^{Qi)/^{Q2) ~ 1-4), and thus it is not needed to 
take care of the anisotropy for a FSS analysis in our spin 
model. 

Figure 7(b) shows the system size dependence of 
the correlation length. To improve the statistical accu- 
racy, we use the average of three correlation lengths 
(C = iUli) + ^(92) + S.ils))/^)- The every correlation 
length shows rapid growth at the transition temperature, 
and they exceed system sizes in the low-temperature 
phase. In this study, every simulation shows markedly 
large correlation lengths below the transition temper- 
ature, and thus we cannot examine whether the true 
long-range order exists or not. However, considering the 
Mermin- Wagner theorem, we expect that the correla- 
tion lengths are finite at finite temperatures even though 
they are extremely large in the low-tcmpcraturc phase. 

3.2 Spontaneous Breaking of Threefold Symmetry 

From the ground-state properties of the J1-J3 model, 
we expect that the spontaneous breaking of the three- 
fold symmetry occurs at finite temperatures. At the 
ground state, along one of the three axes, the inner 
product of NN spin pairs is a negative value charac- 
terized by k (cos(57r/9) ^ -0.17364818), while along 
the other axes, they have a positive value given by k/2 
(cos(57r/18) ^ 0.64278761). Thus, we calculate de- 
fined by 

1 



£u = 



i2 



E 



Si ' Sj 



(/x=l,2,3), (3.8) 



(VJ>NN II axis n 



where index ji indicates one of the three axes of a tri- 
angular lattice (i.e., axis 1, 2, or 3 in Fig. 1). In the 
paramagnetic phase, the values of do not depend on 
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Fig. 3. Temperature dependence of (a) internal energy per site 
E and (b) specific ficat C for lattice sizes L = 36, 72, and 108. 
Insets show enlarged view near the transition temperature. 



/X, because the symmetry does not break. However, in 
the spiral state, one of the e^'s is different from the oth- 
ers. We sort the three averages in descending order and 
define Ei, E2, and as 



El = (max{ei,e2,£3}), 
E2 = (mid{£i,e2,e3}), 
E3 = (min{£i,e2,e3}), 



(3.9) 



where the function "mid{a, b, c}" chooses the second 
largest value from a, b, and c. Since these quantities 
are the bond energies for the NN pairs along each axis, 
we call them the direction-specified bond energies. Fig- 
ure. 8(a) shows the temperature dependence of Ei, E2, 
and E3 for L = 72. In the paramagnetic phase above 
Tc, the direction-specified bond energies take the same 
value, as expected. While below Tc, Ex and E2 increase 
but i?3 decreases. This result implies that the threefold 
symmetry is broken spontaneously at T^. To study this 
anomalous behavior quantitatively, we estimate the en- 
ergy difference defined by Ai? = Ei — E3. The tem- 
perature dependence of AE is shown in Fig. 8(b). The 
energy difference abruptly increases at Tc, and the gradi- 
ent of AE appears to diverge as the lattice size increases. 
Considering the results, we conclude that the first-order 
phase transition accompanies the spontaneous breaking 
of the threefold lattice rotation symmetry. 

Owing to the discrete symmetry breaking, as in the 
Potts model, a finite temperature phase transi- 
tion can occm in the J1-./3 model without violating the 
Mermin- Wagner theorem. Recently, the occurrence of a 
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Fig. 4. Energy distribution P(£;) for L=36(T/J3=0. 4841), L=72 
(r/J3=0.4773), L=108 (T/J3=0.4758), L=144 (T/J3=0.4753), 
and L=180 (T/J3=0.4750). We have confirmed that paramag- 
netic states and spiral states appear with almost the same prob- 
ability in eight independent runs. 
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Fig. 5. First-order phase transition temperature plotted as a 
function of the inverse-square of lattice size. Straight lines in- 
dicate extrapolation of transition temperature. 
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Fig. 7. (a) Temperature dependence of the correlation lengths 
€('72)1 ^^'^ €('73) for L = 72. (b) Temperature dependence 
of the average of three correlation lengths = (^{qi) + £.{<l2) + 
^(q3))/3). Inset shows enlarged view near the transition temper- 
ature. 



first-order phase transition with the threefold symme- 
try breaking has been reported in other two-dimensional 
frustrated continuous spin systems. ^^'^^^ 

3.3 Relation with Z2 Vortex Dissociation 

Since the order parameter space of 120-degree struc- 
ture is isomorphic to P3, the point defect in the AF Ji 
model is ni(P3)=Z2.^-' In the real space, this point de- 
fect corresponds to a vortex configuration labeled by 
or 1, which is called a Z2 vortex. It has been reported 
that the topological transition driven by the dissociation 
of Z2 vortex pairs occurs in the AF Ji model.^'^^'^^^ The 
Z2 vortex configuration is created by two unit vectors, 
(i) One is the vector chirality defined by 

2 

n{r) = — -j= (si X S2 + S2 X S3 + S3 X Si) , (3.10) 
Sv 3 

where the subscript of spin denotes one of the corners 
on an elementary upward triangle at the position r. (ii) 
The pointing vector of one spin is located on an elemen- 
tary upward triangle. In the simulation, the Z2 vortex is 
found by the following procedure. We calculate the rota- 
tion axis n and rotation angle lo between two orthogonal 
coordinates defined by two unit vectors. A variable for 
the link between two elementary upward triangles labeled 
j is defined by 



Uj = exp ■ o") 



(3.11) 



Fig. 6. AF = ln(Pmax(£;)/fmin{-E^)) plotted as a function of lat- 
tice size. Lines between points are visual guides. 
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Fig. 8. Temperature depciideiice of (a) direction-specified bond 
energies Ei, E2, and E-^ for L = 72 and (b) energy difference 
AE (= E1 — E3) for various sizes, fnset sfiows enlarged view near 
transition temperature. 



where <t is the PauH matrix. The vorticity V[c] within a 
closed contour c is defined by 



V[c 



1 




(3.12) 



If the contour c involves a Z2 vortex, V[c] = 1; otherwise 
V[c] = 0. Thus, the number density of Z2 vortices is 
given by 

c 

where the sum is taken over all of the smallest closed con- 
tours c, and Nf. is the number of smallest contours. This 
smallest contour is a x triangle which is formed 
by connections between next NN elementary upward tri- 
angles on the original triangular lattice. 

The order parameter space of the spiral ground state of 
the J1-J3 model is isomorphic to P3, and a Z2 point de- 
fect exists. Thus, the topological transition may occur in 
the J1-J3 model, as well as in the AF Ji model. Accord- 
ingly, there is a possibility of a two-step transition: one 
is the topological transition caused by a Z2 point defect, 
another is the first-order phase transition. To examine 
this possibility, we calculate the number density of Z2 
vortices riy. Since the third NN interaction is dominant 
in the present model, we set twice the lattice constant as 
the unit of length. Thus, the limit of Ji = in the present 
model is equivalent to the AF Ji model. Figure 9 shows 
the temperature dependence of at J1/J3 = r^T^/g. At 
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Fig. 9. Temperature dependence of the number density of Z2 vor- 
tices riv Inset shows enlarged view near the transition tempera- 
ture. 



the transition temperature, the gradient of Uy appears 
to diverge as the lattice size increases. Below the disso- 
ciation temperature, the number density of Z2 vortices 
obeys the Arrhenius law: 



(3.14) 



where 2/x is the chemical potential of the Z2 vortex pair. 
This is because the Z2 vortices are in a bound state be- 
low the dissociation temperature. We show the Arrhe- 
nius plot of Uy for L = 72 in Fig. 10 and find that the 
temperature dependence of riy fits well to eq. (3.14) in 
the low-temperature phase. We estimate the dissociation 
temperature where the data start to violate eq. (3.14). 
At J1/J3 = Tstt/q and —1/3, the dissociation tempera- 
ture coincides with the first-order phase transition point 
within the accuracy of the simulation, and thus we con- 
clude that the dissociation of Z2 vortex pairs occurs at 
the first-order phase transition temperature. 

4. Effects of Easy-Plane Anisotropy 

4-1 Model and Ground-State Properties 

In antiferromagnetic systems on a triangular lattice, 
introduction of easy-plane anisotropy might bring about 
a change in the behavior of the system. For example, 
it has been observed that easy-plane anisotropy causes 
the two separate transitions: the magnetic Kosterlitz- 
Thouless transition and the topological Z2 chirality tran- 
sition.^-'^' ^^'^"'^ Therefore, it is possible that easy-plane 
anisotropy changes the behavior of the J1-J3 model. 

The model Hamiltonian is given by 

n = Ji ^ Si- Sj + Js Si■Sj+D,J2i■^^f^ 



(»J>NN 



(».i>3 



(4.1) 



where sf is the z component of spin Sj and £'z(> 0) is 
easy-plane anisotropy. We study the ground-state prop- 
erties for — 4 < J1/J3 < 0. When easy- plane anisotropy 
is absent {D^ = 0), there are three types of struc- 
tures at the ground state. Since the ground-state spin 
configuration is planar, it is not changed by easy-plane 
anisotropy. However, owing to the term, all the spins 
lie in the XY plane at the ground state, and the order 
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- JylJ^=Q \ \ 

o_ □ JilJo,=-m \ \ 

T /i/y3=-0.73425685 \ \ " 

1 2 3 4 5 
/3/r 

Fig. 10. Arrhenius plot of the number density of Z2 vortices 
for various interaction ratios when the lattice size is L = 72. Lines 
are the results of the least-squares fitting. The specific-heat peak 
corresponding to the first-order phase transition point exists at 
Jg/Tc ^ 2.11 (J1/J3 = rg^/g) and Jg/Tc ^ 2.38 (J1/J3 = 
— 1/3).^"^ The topological transition temperature in the AF Ji 
model is J3/T„ ~ 3 (J1/J3 = o).S. 28,44, 50,51) 

parameter space is isomorphic to the one-dimensional 
sphere Si, which corresponds to U(l) symmetry. The 
spin configurations characterized by wave vector k and 
—k have different topologies in contrast to the isotropic 
case. This is because the two spin configurations are dif- 
ferent from each other under global U(l) spin rotation. 
Thus, in the anisotropic case, there are six distinct groups 
of states corresponding to fe = ±(fc,0), ±(fc/2, -\/3fc/2), 
and ±(fc/2, — •\/3fc/2) at the ground state. The num- 
ber of ground states increases in comparison with the 
isotropic case because of the appearance of topological 
Z2 chirality. Furthermore, the point defect is also dif- 
ferent. The point defect of the J1-J3 model with easy- 
plane anisotropy is ni(Si)=Z, whereas the point defect 
is ni(P3) — Z2 in the isotropic model. In the following 
subsections, we investigate the lattice rotation symmetry 
breaking, the Z2 symmetry breaking of chirality, and the 
dissociation of Z vortex pairs in the anisotropic model. 

^.2 XY Spin Limit 

We investigate finite temperature properties of the XY 
spin limit (Dz — > 00). The interaction ratio J1/J3 is 
set to ^5^/9, which is the same value as we adopted in 
§3.1. Whereas we adopted the heat-bath method for the 
isotropic J1-J3 model, we adopt the Metropolis update 
for the anisotropic model. Figure 11(a) shows the tem- 
perature dependence of the specific heat C. The specific 
heat exhibits a divergent single peak analogous to the 
isotropic case. Figure 12 shows the probability distribu- 
tion of the internal energy P{E) near the temperature at 
which the peak of C is located. Although P{E)'s for each 
system size show a bimodal distribution, the low-energy 
state seldom undergoes a transition to the high-energy 
state and vice versa when the lattice size is 108. Thus, 
we omit the data for L = 108. Since the valley in the 
middle of P{E) oi L = 12 is obviously deeper than that 



of L = 36, we conclude that the phase transition is of 
the first order. 

To investigate the spontaneous symmetry breaking at 
the transition temperature, we calculate the energy dif- 
ference of the direction specified bond energies /S.E and 
the chirality k,. The chirality defined by eq. (3.10) be- 
comes scalar for anisotropic spins. A_E and n are plot- 
ted versus temperature in Figs. 11(b) and 11(c), respec- 
tively. These quantities abruptly increase at the first- 
order phase transition temperature. Ai? and k, corre- 
spond to the order parameters of lattice rotation sym- 
metry breaking and Z2 symmetry breaking of chirality, 
respectively. Therefore, we conclude that the first-order 
phase transition accompanies the spontaneous breaking 
of the sixfold symmetry. 

To study the dissociation of Z vortex pairs, we calcu- 
late the number density of Z vortices n^^ . The number 
density of Z vortices is calculated in the same manner as 
in ref. 52. The Arrhenius plot of n^^ for L = 72 is shown 
in Fig. 13. The data are well fitted by the Arrhenius law 
in the low-temperature phase. We estimate the dissoci- 
ation temperature as described in the previous section. 
In Fig. 13, the dissociation temperature appears to co- 
incide with the first-order phase transition temperature. 
The Arrhenius plot indicates that Z vortex pairs dissoci- 
ate at the first-order phase transition temperature within 
the accuracy of the simulation. From these results, the 
first-order phase transition in the XY spin limit of the 
J1-J3 model accompanies the lattice rotation symmetry 
breaking, the Z2 symmetry breaking of chirality, and the 
dissociation of Z vortex pairs. 

4-.3 Finite Easy-Plane Anisotropy 

In this subsection, we consider a finite easy-plane 
anisotropy {Dz > 0). Figure 14 shows the transition tem- 
perature Tc{L) and the latent heat I as a function of lat- 
tice size L. The transition temperature Tc{L) is defined 
as the point where the bottom of the valley in P{E) di- 
vides P{E) equally. The latent heat / is estimated from 
the width of the bimodal distribution. Both Tc{L) and 
I continuously increase as the anisotropy increases. 
In other words, although the introduction of Dz raises 
the transition temperature and the latent heat, it does 
not change the essential properties of the J1-J3 model. 
Therefore, the Z vortex dissociation is continuously con- 
nected to the Z2 vortex dissociation as the magnitude of 
anisotropy decreases in the J1-J3 model. 

5. Dependence of Interaction Ratio 

So far, we have considered the first-order phase tran- 
sitions of the J1-J3 model. However, when we set Ji 
to be zero, the system does not show the specific-heat 
anomaly.*-* Thus, there should be a threshold value of Ji 
where the phase transition changes from being continu- 
ous to being discontinuous. In this section, we examine 
the interaction ratio dependence of the isotropic Heisen- 
berg model. Since the period of the ground-state spin 
configuration becomes quite large for the small Ji sys- 
tem, we restrict the range of J1/J3 to —1.0 < J1/J3 < 
—0.2. To examine whether the phase transition is con- 
tinuous, we calculate the probability distribution of the 
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Fig. 13. Arrhenius plot of the number density of Z vortices of 
the XY spin limit for L = 72. Lines are the results of least-squares 
fitting. The specific-heat peak corresponding to the first-order 
phase transition point exists at Js/Tc = 1.33 (J1/J3 = ''STr/g)- 
The dissociation temperature of Z vortex pairs in the XY spin 
version of the AF Ji model is J-j/Ty ~ 2 {J1/J3 = 0).''' 
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Fig. 11. Temperature dependence of (a) specific heat C, (b) 
AE (= El — E3), and (c) chirality k of the XY spin limit for 
L = 36 and 72. 
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Fig. 14. Easy-plane anisotropy dependence of (a) first-order 
phase transition temperature Tc{L) and (b) latent heat I. The 
point where = corresponds to the isotropic model. The 
dashed lines indicate the results for the XY spin limit {Dz — > 00). 



internal energy. In the range of J\l Jz used here, the prob- 
Fig. 12. Energy distribution P(£;) of the XY spin limit for L=36 ability distribution is always bimodal near the critical 
(T/ J3=0.7560) and L=72 (T/ J3=0.7498). p^j^^^ ^j^^^^ ^j^^ pj^^^^ transitions are of the first order 

for —1.0 < J\l Jz < —0.2. Figure 15 shows the first-order 
phase transition temperature Tc{L) and the latent heat 
I as functions of J1/J3. The phase diagrams show that 
the transition temperature continuously increases as the 
interaction ratio | J1/J3I increases, and the latent heat is 
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Fig. 15. (a) First-order phase transition temperature Tc{L) and 
(b) latent heat I as a function of — J1/J3 for lattice sizes L = 36, 
48, 60, 72, 84, and 96. The cross points at J1/J3 = indicate 
the temperature of the dissociation of Z2 vortex pairs T^/J:j ~ 
Q 38,28,44,50,51) g^j^jj ^jjg continuous transition (I = 0). 



also proportional to — J1/J3. In the simulation, we ob- 
serve neither a new phase boundary nor a threshold value 
of Ji where the phase transition changes from being con- 
tinuous to being discontinuous. Thus, at least for J1/J3 
to — 1.0 < Ji/ J3 < —0.2, the phase boundary exists only 
between the paramagnetic phase and the spiral phase, 
and the phase transition is of the first order. However, 
we cannot exclude the possibility that the phase transi- 
tion is continuous below a small Ji. We summarize the 
phase diagram of temperature versus interaction ratio 
and easy-plane anisotropy in Fig. 16. Figure 16 shows 
that the first-order phase transition is quite common for 
various environments in the J1-J3 model. 

6. Discussion and Summary 

We discuss the implications of the present results in 
regard to the experimental results for NiGa2S4. In the 
J1-J3 model, even if the value of J1/J3 changes and 
easy-plane anisotropy is added, the thermal phase tran- 
sition is of the first order. This is inconsistent with ex- 
perimental results for NiGa2S4, because clear evidence 
of a first-order phase transition has not been observed. 
Furthermore, the spin correlation estimated from neu- 
tron diffraction measurements reached only a few sites 
at low temperatures, in contrast with the large correla- 
tion length estimated by simulation. Therefore, the low- 
temperature phase of the J1-J3 model are qualitatively 
different from that observed in NiGa2S4. 

In this paper, we have studied the critical phenom- 
ena of the J1-J3 model on a triangular lattice by a 




Fig. 16. (Color online) Schematic phase diagram of temperature 
T versus interaction ratio J1/J3 and easy- plane anisotropy Dz. 
Bold lines indicate the boundaries of first-order phase transition 
found in this work. Blue (light-colored) plane indicates the pre- 
dicted boundary of the first-order phase transition. 



Monte Carlo method. We have set the nearest-neighbor 
interaction Ji and the third nearest-neighbor interac- 
tion J3 to be ferromagnetic and antiferromagnetic, re- 
spectively. When the interaction ratio Ji / J3 is set to be 
—4 < J1/J3 < 0, the spin configuration of the ground 
state forms a spiral. In this spiral configuration, global 
spin rotation cannot compensate for the effect of 120- 
degree lattice rotation, and there are three types of dis- 
tinct structures. We have considered the phase transi- 
tion of the J1-J3 model under three conditions, (i) The 
interaction ratio J1/J3 is set to —0.73425685 such that 
the ground state is commensurate. Adopted spins are 
isotropic, (ii) An extra anisotropic term, -Dz^j(sf)^, is 
added to the Hamiltonian. (iii) Several values of the in- 
teraction ratio J1/J3 are incorporated into the isotropic 
model. 

In §3, under condition (i), we have investigated 
whether the incommensurate nature is essential. Since 
a clear first-order phase transition is seen in the system 
where the ground state is commensurate, we conclude 
that incommensurability is irrelevant to the discontinu- 
ous phase transition. We have calculated the energy dif- 
ference AE and the number density of Z2 vortices riy 
near the transition point. The energy difference AE is an 
indicator of threefold symmetry breaking, and the num- 
ber density of Z2 vortices ny follows the Arrhenius law 
below the dissociation temperature. Since both AE and 
Uy exhibit the anomalous behavior at the same time, the 
breaking of threefold lattice rotation symmetry and the 
dissociation of Z2 vortex pairs occur at the same first- 
order phase transition point. 

In §4, we have considered the effects of easy-plane 
anisotropy on the first-order phase transition under con- 
dition (ii). In this case, we have also found the existence 
of the first-order phase transition, as in the isotropic 
case. Even when easy-plane anisotropy is introduced into 
the model, the order of the thermal transition does not 
change. Furthermore, this first-order phase transition ac- 
companies the lattice rotation symmetry breaking, the 
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Zi symmetry breaking of chirality, and the dissociation 
of Z vortex pairs. Therefore, the essential properties of 
the J1-J3 model are not change by taken easy-plane 
anisotropy into account. In other words, the Z vortex 
dissociation is continuously connected to the Z^ vortex 
dissociation in the Ji- J3 model. 

In §5, under the condition (iii), we have examined 
whether lowering the interaction ratio would bring about 
a continuous phase transition. We have found that a first- 
order phase transition occurs for —1.0 < J1/J3 < —0.2 
in the isotropic case. However, the possibility of the con- 
tinuous phase transition still remains for — 0.2 < J1/J3. 

In spite of the fact that we have adopted several con- 
ditions, as mentioned above, the J1-J3 model always 
exhibits a first-order phase transition. Prom this find- 
ing, the first-order phase transition in the J1-J3 model 
is quite common and robust against changing environ- 
ments. Therefore, wc expect that such a first-order phase 
transition will be observed in actual frustrated magnets 
with competing interactions. 
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